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^ Abstract 

^ The Camassa-Holm equation (CH) is a well known integrable equation describing 

O the velocity dynamics of shallow water waves. This equation exhibits spontaneous 

^ emergence of singular solutions (peakons) from smooth initial conditions. The CH 

equation has been recently extended to a two-component integrable system (CH2), 
p I which includes both velocity and density variables in the dynamics. Although pos- 

sessing peakon solutions in the velocity, the CH2 equation does not admit singular 
^ solutions in the density profile. We modify the CH2 system to allow dependence 

^ on average density as well as pointwise density. The modified CH2 system (MCH2) 

^ does admit peakon solutions in velocity and average density. We analytically identify 

the steepening mechanism that allows the singular solutions to emerge from smooth 
^ spatially-confined initial data. Numerical results for MCH2 are given and compared 

0\ with the pure CH2 case. These numerics show that the modification in MCH2 to 

^ introduce average density has little short-time effect on the emergent dynamical 

Q properties. However, an analytical and numerical study of pairwise peakon inter- 

j>I actions for MCH2 shows a new as^^mptotic feature. Namely, besides the expected 

soliton scattering behavior seen in overtaking and head-on peakon collisions, MCH2 
also allows the phase shift of the peakon collision to diverge in certain parameter 
regimes. 
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1 Introduction 

Singular, mcasurc-valucd solutions in fluids appear in the familiar example of the point 
vortex solutions for the Euler vorticity equation on the plane. The point vortices are 
measure-valued solutions whose motion is governed by a type of multi-particle dynamics. 
In three dimensions this concept extends to vortex filaments or vortex sheets, for which 
the vorticity is supported on a lower dimensional submanifold (ID or 2D respectively) of 
the Euclidean space M^. These singular vortex solutions form an invariant manifold. That 
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is, they persist under the evolution of the Euler inviscid fluid equations if they are present 
initiaUy, but they are not created by Euler fluid motion from smooth initial conditions. 
Whether they may be created from smooth initial conditions by Navier-Stokes viscous 
fluid motion remains a famous open problem. Invariant manifolds of singular solutions 
also exist in plasma physics as magnetic held lines in magnetohydrodynamics and in 
Vlasov kinetic theory as single-particle trajectories. 

Shallow water theory introduces in the limit of vanishing hnear dispersion another class 
of singular solutions which has the property of emerging from smooth initial conditions. 
This class of emergent singular solutions is the main subject of the present paper. 



1.1 Singular solutions for unidirectional shallow water waves 

The Korteweg-de Vries equation (KdV) for unidirectional shallow water waves is 

Ut + 3uUx = -CqUx + 'yUxxx , (1-1) 

in which the fluid velocity m is a function of time t and position x on the real line. 
The subscripts denote the corresponding partial derivatives, while the constants Cq and 7 
represent the effects of hnear dispersion. KdV appears at hnear order in an asymptotic 
expansion for unidirectional shallow water waves on a free surface under gravity. The 
expansion is made in terms of two small dimensionless ratios for small-amplitude long 
waves in shallow water. 

At quadratic order in the same asymptotic expansion, the Camassa-Holm equation 
(CH) appears [3j, 

mt + um^ + 2mu:r = -CqUx + '-^u^xx , m = u- a^u^x- (1-2) 

In the limit that 0, CH recovers KdV. 

Both KdV and CH are completely integrable bi-Hamiltonian equations describing the 
effects of nonlinearity and dispersion on unidirectional shallow water waves at their respec- 
tive orders in asymptotic expansion. Moreover, both KdV and CH arise as compatibility 
conditions for their respective isospectral eigenvalue problem and a linear evolution 
equation for the corresponding eigenfunctions. The properties of being bi-Hamiltonian 
and possessing an associated isospectral problem imply that the one-dimensional KdV and 
CH equations are each completely integrable as Hamiltonian systems. In particular, 
they each possess an infinity of conservation laws and each is solvable by its corresponding 
inverse scattering transform (1ST). Perhaps not unexpectedly, the isospectral eigen- 
value problem for CH differs from that for KdV but recovers the KdV case for — > 0. 



In the absence of linear dispersion the CH equation (1.2) is given by 

mt + urux + 2mux = where m = u — a'^Uxx , (1-3) 

This limit of the CH equation retains only the effects of nonlinear dispersion on unidirec- 
tional shallow water waves and it is still completely integrable. The soliton solution for 
CH without linear dispersion is the peakon defined by 

M(x,t) = ce-l"-'^*!/", (1.4) 



Figure 1: Under the evolution of the dispersionless CH equation (1.3 1, an ordered wave train of 
peakons emerges from a smooth locaHzed initial condition (a Gaussian). The spatial profiles at successive 
times are offset in the vertical to show the creation and evolution of the peakons. The peakon wave train 
eventually wraps around the periodic domain, thereby allowing the leading peakons to overtake the 
slower peakons from behind in collisions that conserve momentum and preserve the peakon shape but 
cause phase shifts in the positions of the peaks, as discussed in j3]. This soliton behavior is the hallmark 
of completely integrable Hamiltonian partial differential equations. The heights/speeds of the solitions 
are governed by the eigenvalues of the associated isospectral eigenvalue problem. 



which is the sohtary travehng wave solution for (1.3). The peakon traveling wave moves 



at speed c equal to its maximum height, at which it has a sharp peak (jump in derivative). 
The spatial velocity profile e"'^'/" is the Green's function for the Helmholtz operator 
(1 — a^d"^) on the real fine with vanishing boundary conditions at spatial infinity. In 
particular, it satisfies 

(1 - a292)e-l^-"*l/" = 2aSix - ct) , (1.5) 



in which the delta function 6{x — q) is defined by 



/(?) = / fix)S{x-q)dx. 



(1.6) 



for an arbitrary smooth function /. 



In its dispersionless limit in (1.3), the CH equation admits solutions representing a 
wave train of peakons 



ulx 



N 



(1.7) 



a=l 



Such a sum is an exact solution of the dispersionless CH equation (1.3) provided the 
time-dependent parameters {pa} and {qa}, a = 1, . . . , N, satisfy certain canonical Hamil- 



tonian equations that will be discussed later. By equation (1.5), the peakon wave train 



(1.7) corresponds to a sum over delta functions representing the singular solution in 
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momentum, 

N 

m{x,t) =2aJ2Pait)S{x-qa{t)). (1.8) 

a=l 

A remarkable feature of these singular solutions is that they emerge from any spatially 
confined initial condition, as shown in the sequential plots of velocity profiles in Figure 
[TJ Being solitons with no internal degrees of freedom, the singular solutions interact by 
scattering elastically with each other. These elastic-collision solution properties hold for 
any Green's function or convolution kernel K{x) in the convolution relation u = K * m 
between velocity u and momentum m [8]. The generalization 

rrit + um^ + 2mUx = where u = K * m , (1-9) 

for an arbitrary Green's function K{x) is called EPDiff, because it arises as an Euler- 
Poincare (EP) equation from Hamilton's principle for a Lagrangian defined as a metric 
on the tangent space of the diffeomorphisms (Diff) [12]. In particular, EPDiff describes 
geodesic motion on the diffeomorphisms with respect to the metric associated with the 
Green's function K. 

The EPDiff formulation generalizes to higher spatial dimensions immediately. In 
higher dimensions EPDiff arises in applications such as turbulence where it is the ba- 
sis for the Navier-Stokes-alpha model [7] and in imaging [HI [TT] where it appears in the 
optimal control approach to template matching. In any number of spatial dimensions, 
EPDiff admits the spontaneous emergence of singular solutions from confined initial con- 
figurations. In ID, this emergent singular behavior results from the steepening lemma, 
proved for dispersionless CH with K{x) = e~\^\l^ in [3]. 



1.2 Plan of the paper 

The dispersionless Camassa-Holm equation has a two-component integrable extension 
(CH2) [H [TBI E]- The CH2 system of equations involves both fiuid density and momen- 
tum, but it possesses singular solutions only in the latter variable. Section [2] discusses 
CH2 in the context of shallow water systems. After discussing some of the shallow-water 



properties of CH2 in Section |2.1[ a modified system MCH2 is proposed in Section [2^ 
Although the MCH2 system may not be integrable, it does allow delta-like singular solu- 
tions in hoih variables, not just the fiuid momentum. These singular solutions exist in any 
number of dimensions, as discussed in Section |3| The mechanism of nonlinear steepening 
by which these singular solutions form is explained via a lemma in Section |4} This lemma 
reveals the conditions under which wave breaking occurs, so that the MCH2 fiuid velocity 
in one spatial dimension develops a negative vertical slope, starting from smooth initial 
conditions. Section [5] presents numerical results of emergence and stability of MCH2 solu- 
tions for a typical fiuid problem known as dam breaking. The numerical results indicate 
that the singular solutions are quite stable as they emerge and when they collide, but may 
tend toward the pure CH peakon solutions asymptotically in time. Section [6] analyzes the 
pairwise collisions of MCH2 peakons in one dimension. Section [7] addresses the numerical 
phenomenology of attractive interactions of peakons, in which the sign of the gravitational 
acceleration is opposite that of the shallow-water case. Section [8] concludes with some 
remarks about the potential applications of MCH2 in imaging science. 
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2 Two- component systems 

2.1 The two-component Camassa-Holm system (CH2) 



In recent years, the Camassa-Holm (CH) equation (L3| has been extended so as to com 



bine its integrability property with compressibility, or free-surface elevation dynamics in 
its shallow-water interpretation. This extension involves adding a continuity equation for 
the scalar density (or total depth) p G for real functions T and including a pressure 
term involving p in the equation for the fluid momentum, as weU as the fluid velocity u. 
The CH2 equations in that case are specifically [H [6l [181 123] 



rat + ura^ + 2mUx = —g ppx where m = u — o?u 



XX 1 



Pt + {pu)x = 0, (2.1) 

where g > is the downward constant acceleration of gravity in applications to shallow 
water waves. Boundary conditions are taken as u and p — > po = const as |x| oo. 



Complete integrability of the CH2 system (2.1) may be proven by writing it as the 



compatibility condition for two linear systems with a spectral parameter A as in [H [TSl E] 

ipxx = + — (7p^A^j?/' (Eigenvalue problem) (2.2) 

ipt = (wT^u^i^x + T^Uxip (Evolution equation) (2.3) 
V 2A / 2 

Requiring compatibflity ipxxt = "iptxx and isospectrality dX/dt = recovers the CH2 system 



(12T|. 

Geometrically, the CH2 system corresponds to geodesic motion with respect to the 
conserved metric [HI [IT] 

L{u,p) = ^||n||^i + |||p-po||i2 (2.4) 

for velocity and density {u, p) G TDifi^ x JF, where JF denotes the space of scalar functions. 
The cases g > (resp. g < 0) correspond to repulsive (resp. attractive) interactions 
among particles of positive pointwise density p. The CH2 system emerges from Hamilton's 
principle 6S = with the action S = J L{u, p)dt when the Lagrangian L{u, p) is taken to 
be the metric in (2.4). The CH2 system ( 2T| has been shown to possess peakon solutions 



in velocity u and corner-like solutions in density p in |4j. However, singular solutions 
do not exist for its density variable [5j. Singular solutions in density p (and therefore 
reduction to a finite-dimensional system of Hamiltonian equations) wfll be restored in the 
next section by a shght modification of CH2. 

2.2 Modified CH2 (MCH2) 



The proposed modification of the CH2 system (2.1) is expressed in terms of an averaged 



or filtered density g in analogy to the relation between momentum and velocity by setting 



m = (1 — ald'^) u and p = (1 — a^d'^) g , 
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where one denotes d = d/dx and defines two length scales cti and 02. This modification 
will amount to strengthening the norm for 'q from to in the potential energy term 



in the metric Lagrangian (2.4) in Hamilton's principle for the CH2 system. 



The modified CH2 system (MCH2) is written in terms of velocity u and locally- 
averaged density 'q (or depth, in the shallow- water interpretation). MCH2 is defined 
as geodesic motion on the semidirect product Lie group Diff (§) with respect to a certain 
metric and is given as a set of Euler-Poincare equations on the dual of the correspond- 
ing Lie algebra X(S)JF. In the general case, for a Lagrangian L{u,g), the corresponding 
semidirect-product Euler-Poincare equations are written as [12] 



(2.5) 



dt 6u ^ 6u 6g ' dt Sg ^ 6g ^ 

where £u (SL/Su) is the Lie derivative of the one-form density m = 6L/Su with respect to 
the vector field u and £uSL/6g is the corresponding Lie derivative of the scalar density 
SL/Sg. 



The integrable CH2 system and thee modified system MCH2 may both be derived as 



semidirect-product Euler-Poincare equations (2.5) from the following type of variational 
principle defined on the Lie algebra X(§)J-' 

S [ L{u,g)dt = 0, 

J to 



with Lagrangian 

L{u,g) 



\u 



— I|2 

QoWm 



2 2 2 \ 1 
u + aiU^) ax 
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dx . 



(2.6) 
(2.7) 



in which the last fine defines the norms of m and (g — go) and gQ is taken to be constant. 
The variational derivatives of this Lagrangian define the variables m and p as 

9{l-ald'){-g ^ 



5g 
5L 
5u 



(1 — a^d )u =: m . 



9P: 



Substituting these variational derivatives into the Euler-Poincare equations (2.5) recovers 
the CH2 equations (2.1) for the constant values g > 0, ai = a and 0^2 = 0. However, 
when a| > the MCH2 equations result. Namely, in one dimension the MCH2 system is 



-9PQx^ 
0, 



rrit + UTJix + 2mUx 
Pt + {pu)x 

where m and p are defined in terms of u and 'g via the Helmholtz operators 

m := {1 — a\d'^)u and p := {1 — a\d'^){~g ~' 



f2.8) 



(2.9) 



This shght modification of CH2 defines the MCH2 equations in one dimension and, as we 
shall show, it suffices for the existence of emergent singular solutions in both m and p. 
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Remark 2,1 (CH2 and MCH2 in higher dimensions) In higher dimensions the La- 
grangian (2.1] generalizes to 

L{nrQ) = \ Ju- (l-a?A)ud"x+| J {g -g,){l - a'.A) {g -g,) d"^^ , (2.10) 

where A is the n- dimensional Laplacian operator. The semidirect- product Euler-Poincare 
equations for this Lagrangian produce the desired generalizations of CH2 and MCH2 to 
higher dimensions. 



3 Existence of singular solutions 



One evaluates the variational derivatives of the n-dimensional Lagrangian (|2.10| as 

m ■ dx ® dV^ and 



5L 
5u 



SL 

6g 



gp0dV, 

■jg 

with co-vector m and scalar function p defined by 

'1 — «^A)u and p=(l — a2A)^. 



(3.1) 



m := 



Then the semidirect-product Euler-Poincare equations (2.5) in 
coordinates as 



(3.2) 

°" may be written in 



dm 

~dt 

dp 
di 



Convection 
div pu . 



(Vu) 



m 



Stretching 



m(div u) 
Expansion 



gp'^g 

Force 



(3.3) 



(3.4) 



Here one denotes (Vu) -m = rrijVuK To explain the terms in underbraces, we rewrite 
EPDiff as the change in the one-form density of momentum along the characteristic curves 
of the velocity. In vector coordinates, this is 



d 
dt 



m • c/x ® dV 



gpdg ® dV along 



(ix 

Itt 



u = i^i * m . 



(3.5) 



This form of the EPDiff equation also emphasizes its nonlocality, since the velocity is 
obtained from the momentum density by convolution against the Green's function Ki of 
the first Helmholtz operator in (3.2). Likewise, the average density 'g follows from the 
pointwise density p by convolution against the Green's function K2. 

A direct substitution into the Euler-Poincare equations (3.3) and (3.4) in M" coordi- 
nates shows that these equations possess the following singular solutions 

N 



{in,p)=J2 I (P.(s,t),Wi(s))5(x-Q,(s,t)) dh. 
i=i ^ 



(3.6) 
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where (Qi,Pj) satisfy 



dt 



dt 

Ms)w,{s') ^K,{Q,{s,t)-Q,{s',t)) dV, 

where Ki and K2 are the kernels for the Green's functions of the Helmholtz operators 
with length scales ai and a2, respectively. The continuity equation (3.4) imphes that 
the weights Wi are independent of time. These equations may be written as Hamilton's 
canonical equations with Hamiltonian 

1 ^ rr 

Hn = ^Y. rPr{s,t)-V,{s\t) K,{Cl,{s,t)-Qlj{s\t)) d'^sdV 



2 



N 

+ ^^^^is)w,{s')K,{Q,{s,t)-Q,{s',t)) d'^sdV, 



2 

Remark 3.1 (Why are these equations canonical?) The explanation of why the dy- 
namical equations for the singular solutions are canonically Hamiltonian is beyond the 
scope of the present paper. This is explained from the perspective of Lie group actions and 
momentum maps in fW^ . 

In the one dimensional treatment considered in the remainder of this work, the equa- 
tions above simplify to 



dt 
dF 



dP d d 



(3.7) 

and if we choose a2 = 1 so that K2{x) = ^e"'^', we recover the peakon solution for g 



-g{x,t) = {i-dir'p = lYl 



as well as the usual peakon wave train solution (1.7) for the velocity u. 

Remark 3.2 (Singular potential terms and integrable cases) In the limiting case 
a2 = 0, one recovers the CH2 system corresponding to geodesic motion with a delta- 
function interaction potential. In terms of particle solutions, this means that for a positive 
potential {g > 0), two particles will bounce off immediately before colliding; while for a 
negative potential {g < 0), they will proceed together, attached one to the other and they 
will never split apart. This kind of singular delta-like potential is also present in another 
integrable system, namely the Benney equation /H 
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4 Steepening lemma for MCH2 

Do the singular MCH2 solutions emerge from smooth initial conditions? This is not 
automatic. For instance, we know that point vortices are invariant singular solutions that 
are not generated by vorticity dynamics with smooth initial conditions, whereas peakon 
solutions are always spontaneously generated by the CH flow. 
This section and the next address the foUowing two questions: 

1. Under what conditions does MCH2 admit emergent singular solutions? 

2. What are the stability properties of such solutions? 

This section proves a steepening lemma which reveals the conditions under which the 
singular solutions for MCH2 in one spatial dimension emerge from smooth initial condi- 
tions. 



On using the definitions of m and p in (2.9), the MCH2 equations (2.8) may be 



rewritten in a form that transparently displays its nonlocal nature 

ut + uu^ = -d{Ki*p) with p := u"^ + Y ul + -g"^ 2^^^' ^^'^^ 

Qt + ug^ = -K2* (al {ua:g^)^ + (u^g)^ . (4.2) 
where Ki and K2 are the Green's functions for the corresponding Helmholtz operators. 



For g > 0, conservation of the norm (2.6) ensures that u and g are bounded pointwise. 



because of the Sobolev inequality vahd in one dimension that, 

max (u'^{x,t) + g(g{x,t) — gQ^) < \\u\\'jji + g\\g — gQW^^i = const . (4.3) 

Lemma 4.1 (Steepening Lemma for the MCH2 system) 

Suppose the initial profile of velocity u{0, x) has an inflection point at x = x to the right 
of its maximum, and otherwise it decays to zero in each direction. Under MCH2 dynamics 
with g > a suffi,ciently negative slope at the inflection point will become vertical in finite 
time. 

Proof. Consider the evolution of the slope 

s{t) = u^{x{t),t) , 



at the inflection point x = x{t). The spatial derivative of the w-equation (4.1) yields an 
equation for the evolution of s. Namely, using Uxx(x(t),t) = the spatial derivative leads 
to 



^ + ^2 = -dl{Ki*p) with p 



2 «i 2 , 5^-2 9^1 -2 

u +y3 +-^g -^gx 



x=x(t) 



a\dl)K^*p-\K,*p 



^jp-^Ki*p. (4.4) 



a^ a 
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This calculation implies with g > 



ds 
di 



a 



< 



< 



1 2 , 9(4 

2 "^-^ 



2ai 



2 9^1 



1 



a 



oo 
oo 



-oo 
oo 



-\x-y\/ai 



-\x-y\/ai 



2.9—2 



9^2 -2 
9 



u +-9 - 



9(^2 



x=x{t) 



'Qldy + 



2ai L 



u 



9-2 

- Q 
2 ^ 



x=x{t) 



Qofx dx + 



at 



2.9—2 
U +-0 



(4.5) 



x=x{t) 



where the second step drops the negative terms in the previous hne and the last step uses 
g-k-j/l < X By conservation of the norm with g > in the Sobolev bound (4.3), the 
sum of last three terms remains finite, say less than a number M/2. Consequently, we 
have 

ds 
dt 



< 



1 . M 
-s^ H 

2 2 



which implies, for a sufficiently negative initial slope, s < 
negative and becomes vertical in finite time, as 



(4.6) 

M , that the slope remains 



s{t) < VMcoth 0- 



(4.7) 



where a is a negative constant that determines the initial slope, also negative. Hence, at 
time t = — 2(j/a/M the slope becomes negative and vertical. This wave-breaking result 
for the slope of the fiuid velocity proves the steepening lemma for the MCH2 equation 
with g > 0. ■ 



Remark 4.2 The calculation corresponding to (4.5) for g < results in the estimate 



ds 
dt 



< 



\9\al 
a? 



g'^ dx + 



\9\al 



- x=x{t) 



However, without any control on the pointwise value of g^, no real conclusion may be 
drawn from this estimate. 

Remark 4,3 (Density at the point of wave breaking) One may ask how the den- 
sity behaves when the velocity slope at the inflection point becomes vertical in finite time 
for g > 0. This may be determined by integrating the continuity equation for p along a 
Lagrangian path x = x{xo,t) with x(yXo,0) = xq. Evaluating at the inflection point yields 

1 dp{xo,ty 



p{xo,t) dt 



-sit) 



x{xo,t)=x(t) 



SO that the density of the Lagrangian parcel currently occupying position x(t) is given by 



p{xo,t) 



x{X{),t)=x{ 



= p{xo, 0) exp ^ - ^ 



s{T)dT 



Because the blow up in s{t) as coth in equation (4.I] is integrable, the density of the 
Lagrangian parcel at the inflection point remains finite even though the velocity develops 
a vertical slope. Thus, wave breaking in the fluid velocity does not imply singularity in the 
pointwise density at the point of vertical slope. 
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Remark 4.4 (The formation of delta singularities in momentum) 

Peakons have no inflection points, so they are exempt from the steepening lemma. As 
shown in Figure [I], their formation creates a new inflection point of less negative slope. 
The negative slope of the inflection point present in the initial velocity profile steepens and 
the velocity at the inflection point rises until it lifts above the initial maximum velocity. 
Then it propagates away, leaving an inflection point behind, which will soon repeat the 
process. At each peak in the velocity profile the slope takes a jump, so the curvature there 
is infinite. This in turn leads to the delta singularities in the momentum variable (m) 
obtained by applying the Helmholtz operator to the velocity profile. Finally, the formation 
of the delta-function singularity in momentum m must also be accompanied by a delta- 
function singularity in pointwise density (p), as shown by direct substitution into the 
MCH2 equations in Section\^ Thus, the reaction to wave breaking by nonlinear steepening 
at inflection points of negative slope creates peakons and is the mechanism under which 
MCH2 forms its delta- singularities in m and p. 

5 Numerical results benchmark problem: dam break- 
ing and bores 

The simulations reported in this section solve the following equations with g = 1 

dtu = Ki* {-3uu^ + 2a^M^Ma;a; + ajuu^xx - 9{q - a2Qxx)Qx) > 

dtg = -K2*d{{g-alg^Ju) , (5.1) 

where ai and 02 are the two length scales of nonlocality. The first equation may also be 
written in conservative form as 

ut + uu.j, = -Ki* d \u + —u^ + -Q Y^"") ■ 

Periodic boundary conditions are imposed on the domain [— L, L] with dam-break initial 
conditions given by 

u (x, 0) = 0, g {x, 0) = 1 + tanh(a; + a) — tanh(x — a) , (5.2) 

where a <^ L. 



5.1 Dam breaking for g > Q 

The dam-break arises when a body of water of uniform depth is retained behind a barrier, 
in this case at x = ±a. If this barrier is suddenly removed at t = 0, then the water fiows 
under gravity. The problem is to find the subsequent fiow and moreover, to determine 
the shape of the free surface. This question is addressed in the context of shallow-water 
theory by Acheson (see fTj), and thus serves as a typical hydrodynamic problem to be 
discussed in our framework. 

In the following paragraphs, we compare the numerical results obtained in two differ- 



ent cases: the exactly integrable CH2 system; and the modified MCH2 equations (2.8) 
proposed above (with 02 7^ 0). 
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Case 1: pure CH2 equations. We take ai = 0.3, 02 = 0, a = 4, and L = 127r. The 
choice a2 = gives a density equation (the equation in g) that is unsmoothed. As shown 
in Figure [2} the singular solutions emerge after finite time, and both variables have this 
singular property. More and more singular elements appear as time progresses. 





-20 20 



-20 20 



(a) 



(b) 



Figure 2: (Color onhne) Dam-break results for the pure CH2 system show (a) evolution 
of the density p = Q and (b) evolution of the velocity u from initial conditions (5.2). 



In contrast, one may compare the evolution of the velocity u for the single-component 
CH equation launched from a spatially confined velocity profile: u{x, 0) = tanh(a; + 4) — 
tanh(a; — 4) and g{x, 0) = const with ai = 0.3. In contrast to CH2, only rightward-moving 
peakons emerge from CH as illustrated in Figure [2] Thus, the limit ^2 — > breaks the 
left-right symmetry of the CH2 system. 




-15 -10 -5 5 10 15 



Figure 3: (Color online) Evolution of the velocity u for single-component CH from the 
spatially confined initial velocity profile, u (x, 0) = tanh(x + a) — tanh(x — a). 



Case 2: MCH2. We take ai = 0.3, 02 = 0.1, a = 4, and L = 127r. The solution 
behavior in this case for the modified system with small ^2 shows only shght difference 
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from the pure CH2 case, as seen in Figure [4| One difference is that now the singularities 
have a particle-hke nature, whose collective dynamics can be studied using Hamilton's 



canonical equations (3.7). 




li) 20 




(a) 



(b) 



Figure 4: (Color online) The appearances of the dam-break results for the MCH2 and 
CH2 systems are identical. See Figure [2} The differences plotted in this figure show (a) 
evolution of the filtered density difference q—Qq and (b) evolution of the velocity difference 
u — Uq. Here the subscript zero refers to CH2, for which a2 = 0. 



/ 





(a) 



10 




(b) 



Figure 5: (Color online) Reducing the size of the periodic domain for the dam-break for 
the MCH2 system allows the peakon wave train to wrap around, thereby allowing the 
leading peakons to suffer head-on collisions with the traihng ones. One observes a phase 
shift, during which the peakons form an unstable intermediate state, which then splits 
into peakons again. This phase shift is discussed further in Section [61 



Remark 5.1 One might question whether the Steepening Lemma 4J_ for MCH2 applies 
to the dam-breaking initial conditions, because a zero-velocity initial condition has no 
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inflection points. However, the numerics shows that as the dam breaks the right-hand- 
side of MCH2 immediately generates nonzero velocity. The velocity must develop at least 
one inflection point of negative slope to remain confined with finite speed. This situation 
qualifies it for the lemma. The lemma says that once an inflection point of negative slope 
develops, wave breaking has a deadline. The MCH2 equation avoids wave breaking by 
successively forming peakons, each time shifting the infl,ection point backward before its 
deadline and restarting its clock at the new location. 



5.2 Bore 

An alternative to the stationary dam breaking problem is to consider the dam breaking 
initial condition as moving in a frame of uniform motion and regard the result as modelling 
the development of a bore produced by a moving elevation of water. 



Equations (5.1) are solved with g = I and the following alpha-values: ai = 0.3 and 



«2 = 1, together with the initial conditions 

u (x, 0) = tanh (x + 4) — tanh (x — 4) , ^ (x, 0) = 1. 

The larger value of a2 was needed to prevent gradients in g from becoming nearly vertical 
for the chosen resolution. Gradients in ^ are not a problem for the dam-break initial 
conditions because the traihng density peaks in that case remain small in magnitude. 
However, in the bore configuration the density builds up as time proceeds. Since a large 
value of a2 was needed to control the gradients in g for the bore, a comparison with the 
^2 = case was not possible. 



Remark 5.2 Similar results are obtained for various other initial density profiles. For 
example, with g{x,0) = 6 [tanh (x + 4) — tanh (x — 4)] and b = 0.1, the density profile 
remains finite, and the velocity profile is very similar to that in Figure\^for CH. Similarly, 
for g (x, 0) = 0, the density profile remains zero for all time, and the velocity profile again 
closely resembles that in Figure\^ 



6 Pairwise collisions in one dimension for ^ > 



The interaction of two singular solutions may be analyzed by truncating the sums in (3.6) 
to consider N = 2. In one dimension this yields 

H2 = ^(^P^ + P| + 2 Ki(Qi - Q2) P1P2) + ^-(^wl + wl + 2 K,{Qi - Q2) W1W2) 
and we choose g = 1, corresponding to repulsive interactions. Following [T3], one defines 

P = Pi + P2 ; Q = Qi + Q2] p = P1-P2, g = Qi - Q2 ; W = wi + w2; w = wi 

Consequently, the Hamiltonian can be written as 
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Figure 6: (Color online) (a) Space-time evolution of the function g; (b) Space-time evo- 
lution of the velocity u for the MCH2 system in the bore configuration. Figure [6)d for 
the velocity of the MCH2 system is quite similar to Figure [3] for CH, but close inspection 
shows that it is not identical. 



At this point one writes the canonical equations 
0, 



with W = w = 0. These equations possess the first integral, 

2 



dP 


dn 


~dt " 


dQ 


dp 


--2— 


dt ~ 


dq 



(6.1) 



(^) " " + i^V'' - (1 - /<■,(,)) ] (1 - A'i(,)) 



which finally produces the quadrature 

dt = - 



dKi 



K[^P^ (K^iq) -l)+[An-2W^ + {W^ - w^) (1 - K2{q)) ] (1 - K,{q)) 
Now, choosing Ki = K2 = K and defining the asymptotic speeds ci and C2 such that 



P = Ci + C2 



H= -{cl + cl + wl + wl) 



leads to the momentum relation 

4(ciC2 + W1W2) 



2 I 2 
p + W 



+ {Wl + W2f + {Ci + C2f . 



l-K 

If the smoothing kernel is chosen to satisfy K{Q) = 1 (e.g. inverse Helmholtz or Gaussian), 
then the quantity p^ + w"^ (and thus the momentum p, since w is constant) diverges for 
g = 0. Defining the vector 

= {Ci,Wi) 

allows the problem to be classified into two separate cases: 
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Overtaking collisions: These satisfy rj^ ■ r]2 > 0. Consequently, the peak separation 
cannot vanish, since p"^ + w"^ > 0. In this case it is possible to express the minimum 
separation between the peaks by simply setting = 0: 

r.^ ^ . 4(C1C2 + W1W2) 

(Ci + C2)'' + 4:WlW2 

To visualize such a collision, we take the initial conditions 

qi (0) = 20, (0) = -0.4, wi (0) = = 1, 
q2 (0) = -20, p2 (0) = 0.4, W2 (0) = ^2 = 1, 



and evolve the two-particle equations (6.1) forward in time, using Helmholtz kernels with 
ai = ^2 = 3 (these are the values we use throughout this section). The results are shown 
in Figure [7| The particles initially move towards each other, they interact and feel their 




Figure 7: Particle trajectories for 77^^ ■ > 0. 

mutual repulsion, and then move away from each other indefinitely. The particles are 
always separated from each other by some finite distance and thus they do not change 
their order. 

Head-on collisions: These possess the property that fj^ ■ < 0. Thus the peaks can 
overlap and the momentum p diverges when this happens. For this particular example, 
we choose to study the special case of two completely antisymmetric peakons, so that 
rji = —r]2 and qi = —q2- This is an interesting simplification since it yields P = Q = 
W = and the quadrature formula becomes 

dt = , = , (6.2) 

K'^[m-wHl-K{q))]{l-K{q)) 

where 

AH -w^l- K{q)) = (1 - K{q)) > . 
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Since we are dealing with repulsive interactions, the two particles travel oppositely 
apart after the colhsion. Thus, asymptotically K{q) — 0, p — 2c and H \r]f = cl + e^, 
where rj = (cq, Cq) (or tj = — (cq, Cq)) collects the asymptotic speed and the charge of the 
particle (or antiparticle). Setting if = Cq + Cq in the equation above yields the phase 
trajectory 

SO that p diverges when g — > 0, thereby recovering the same qualitative behavior as the 
ordinary peakon solution for the CH equation. 

We obtain a numerical solution in this case by focussing on the intial conditions 

qi (0) = 20, pi (0) = -0.4, wi (0) = wi = 1 , 

q2 (0) = -20, Co = P2 (0) = 0.4, -cq = W2 (0) = W2 = -I . 

The evolution of these initial conditions is displayed in Figure [8| The particles initially 
move towards each other, collide, repel each other, and separate. Notice however that 




Figure 8: (a) Particle trajectories for rji ■ r]2 < 0. The difference q = qi — q2 changes 
sign at the collision point and, thus, the particle orders are reversed; (b) Phase portrait 
p versus q for this colhsion. 

the particle orders change at the collision point: the particle with label 1 is orignally on 
the right but after the colhsion is on the left, and conversely for particle 2. This is in 
contrast to the case of the overtaking colhsions (case 1). A similar order-changing result 
is observed for the Gaussian kernel. 

Superficially, the idea that the particles exchange order at the collision time may 
appear undesirable, and could be rectified by exchanging the particles' momenta after 
the colhsion. This is legitimate because any choice of method for integrating the equa- 



tions (6.1) beyond the collision singularity is not unique. Nevertheless, the order-changing 
conditions are more natural in this context. To see this, we regularize the kernel K (g) 
by letting K (q) — > i^g (g) = e"*^'^!"'"'^^/", so that (g) is strictly less than unity. In this 
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case, the momenta are finite at the colhsion, and the particle trajectories pass through 
one another. This indicates that the order-changing conditions on the momenta in the 
unregularized case are the more natural choice for integration beyond the collision singu- 
larity. 

We gain some understanding of the mathematical structure of the colhsion by studying 



the quadrature formula (6.2), which in terms of the initial conditions (co,eo) reads 



dt 



dq 



with solution 
where 



^[^{cl + el)-Ael{l-K{qm^-K{q)) 



t = (t){q)~(t) (go) , g (0) = go > 



a 

^0 



a 



Co 



log [cl + el) 



log 



2c^ 



-q/a 



e --c^ + e2e"'?/" + 2co 



The collision time is obtained by setting g = 0, 

tc = la log {cl + el) - 0(go) , 

which is finite. For the starting conditions used in Figure [8| we obtain a collision time 
tc = 47.77, in agreement with the numerical results shown in the figure. Finally, Figure 
[s] (b) provides further information about the collision through a phase portrait; this plot 
coincides exactly with formula (6.3). The momentum is p = — 2co = —0.8 before and 
after the collision, while during the colhsion, it diverges to p = —oo. 



Limiting case: There is also a hmiting case when rj^ ■ r]2 = 0. In this case, one finds 
that assumes the constant value 

p'^ = + -w^ = 4:WiW2 + (ci + C2f . (6.4a) 

Consequently, from the equation 

^^=p{l-K{q)) (6.4b) 

we conclude that if p = — y/ P'^ + — w^, then the two particles merge after the peaks 
overlap and they never split apart. This happens independently of the particular choice of 
K. Other stationary states of the pairwise interaction dynamics may be found by setting 
g^ = and by solving in g. 

To illustrate the dynamics in this hmiting case, we take the following initial conditions: 

gi (0) = 20, pi (0) = -1, wi (0) = ^1 = 1, 



q2 (0) = 



-20, P2 (0) = 1, W2 (0) = ^2 = 1, 
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and evolve the two-particle equations (6.1) forward in time, as shown in Figure |9 



As in the attractive case of head-on collisions, the particles initially move towards each 
other, and interact. In this case, it is easy to understand the mathematical structure of 



the collision because the equations of motion have the simple form exhibited by Eqs. (6.4) 
for which p < 0. Thus, defining r = \p\t, Eq. (6.4b) becomes 



dq 
dr 



Using the indefinite integral 



dq 



we obtain the solution 



1-K{q) 



-a log (1-Coe-^/") 



[l-K{q)]. 



sign(g)alog(l-e-l^l/") 



-90/0 



< Co < 1 . 



The colhsion occurs when g = 0, that is, when r = 00. (It is clear that the collision times 
for the Gaussian and top-hat kernels are also infinite in this case.) Figure |9] thus shows for 
^1 ■ ''72 = that the particles move asymptotically closer to each other, but never actually 
reach the colhsion point. In this case, one may say that the phase shift of the colhsion 
diverges. One might ask how the phase shift behaves in the limit as r]^ - 7)2 ^ 0"*" in a 



sequence of positive values. Figure [10] answers this question by showing that the particles 
approach each other and stay together longer and longer in the hmit as ?7i ■ 0+. 
Figure [To] indicates that the phase shift diverges logarithmically in this hmit. 



60 r 

50- 
40- 



- 30- 




Figure 9: Particle separation trajectories for the hmiting case t)^ ■ = show that the 
particles asymptotically approach each other. 



7 Attractive interactions 



Evolution by either CH2 or MCH2 yields a geodesic flow on the semidirect product group 
Diff(s)jF governed by the Lagrangian (2.7), which is a norm when g > 0, correspond- 
ing to repulsive particles. However, one might also consider, for example, gravitational 
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Figure 10: Particle separation trajectories show a phase shift that lengthens as -772 ^ 
(through positive values). Apparently the phase shift diverges logarithmically as 5 — > 0+. 

interactions, for which particles with mass undergo attractive interactions. Attractive 
interactions correspond to (7 < in the potential term in g of the MCH2 Lagrangian 



(2.7), which is then no longer a norm. It is known that the change in sign of g preserves 
integrability of CH2 when K2 is a delta function [H [H]. One may repeat the previous 
analysis for attractive interactions by considering g < 0. This is the subject of the present 
section. First we show numerical results, then we discuss two-peakon interactions. 

7.1 Numerical results for g < 

To demonstrate the spontanteous emergence of singular solutions for g < 0,we numerically 



solve the equations dealing with attraction. Specifically, we study equations (5.1) with 



g = —1. We solve these equations on a periodic domain [— L, L] with the initial conditions 

(47ra;\ 
— j . 

This kind of initial conditions was also used in [13] to observe "clumpons" in the density 
g. We take ai = 0.3, a2 = 0.1. The results of the numerical simulations are shown in 



Figure 11 These results are very similar to those found in |T3]. However, the types of 



equations involved (Hamiltonian versus gradient fiow) are quite different. 

Remark 7.1 An analytical explanation of the mechanism for the formation of singulari- 



ties shown in Figure \11\ has eluded us, so far. It remains an outstanding problem to prove 
an "anti-maximum principle", that ifg initially has a maximum, then the pointwise p must 
develop a singularity in finite time. Likewise, no theorem exists yet for the formation of 
the velocity singularities. 
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Figure 11: (Color online) Evolution of the system of equations (5.1) for the attractive 
case g = —1. Subfigure (a) shows the smoothed density while (b) plots the velocity u. 
Singularities in the solutions of both variables emerge after finite time. The singularities 
in density p form at the maxima of the initial conditions in g. Positive (resp. negative) 
singularities in velocity u form at the maxima (resp. minima) of the initial conditions in 
g. 

7.2 Pairwise collisions for g < 

The case of attractive potential {K2 = —Ki) yields the following momentum relation 

2 _ 4(ciC2 + W1W2) , ,^^^2 2\ 1 + ^ , / , 



P 



l-K 

Thus, by adding (W'^ — vo^) /{I — K) at each side, one finds 

2 , TI/2 4(ciC2 - W1W2) 2 , ( , n2 

p + W = h w + (ci + C2) . 

i — A 

Consequently, overtaking collisions are defined by C1C2 > W1W2, for which the peaks cannot 
overlap. The head-on collisions take place when C1C2 < WiW2- The momentum formula 
above can be obtained from the case of repulsive interactions by exchanging W2 — > —W2- 
Hence, analysis similar to that in Section |6] yields analogous formulas. Also in the limiting 
case C1C2 = W1W2, one recovers the merging phenomenon already observed in the repulsive 
case. Other stationary states may be found by setting g = 0. In higher dimensions, these 
would correspond to equilibria such as the orbits of binary systems typical of purely 
gravitational interactions {K2 = — A~^). 



8 Applications of CH2 and MCH2 in imaging science 

Much of the discussion here is motivated by the applications of EPDiff in Computational 
Anatomy (CA). The problem for CA is to determine the minimum distance between two 
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images as specified in a certain representation space, V, on which the diffeomorphisms 
act. Metrics are written so that the optimal path in Diff satisfies an evolution equation. 
This equation turns out to be EPDiff', when is a closed contour representing the planar 
shape of the image. The fiow generated by the EPDiff equation transforms one shape 
along a path in the space of smooth invertible maps that takes it optimally into another 
with respect to the chosen norm. Its application to contours in biomedical imaging, for 
example, may be used to quantify growth and measure other changes in shape, such as 
occurs in a beating heart, by providing the optimal transformative mathematical path 
between the two shapes. A discussion of EPDiff and the apphcation of its peakons and 
other singular solutions for matching templates defined by the contours of planar image 
outhnes appears in j.l4] . 




Metamorphosis 

Metamorphosis is a recent development in the problem of image-comparison for CA. In 
the metamorphosis of smooth images one considers a manifold, N which is acted upon 
by a Lie group G: N contains what may be regarded as "deformable objects" and G is 
the group of deformations, which is the group of diffeomorphisms in most applications. 
Several examples for the space N are treated in [T7] . 

Definition 8.1 A metamorphosis is apair of curves {gt, rjt) G GxN parameterized 
by time t, with go = id. Its image is the curve rit & N defined by the action rit = gtTjt- 
The quantities gt and rjt are called the deformation part of the metamorphosis, and 
its template part, respectively. When rjt is constant, the metamorphosis is said to be a 
pure deformation. In the general case, the image is a combination of a deformation 
and template variation. 

Riemannian metric A primary application of the metamorphosis framework is based 
on the definition of a Riemannian metric on G x which is invariant under the action of 
G: {g, ri)h = {gh, h'^rj). The corresponding Lagrangian on TG x TN then takes the form 



where u := gg^^ G Q (the Lie algebra q of the Lie group G), n := 5^77, u := gfj and 
z/)||^ is a norm on g x TN/G parameterized by n. The optimal matching problem 
is now equivalent to the computation of geodesies for the canonical projection of this 
metric from GxN onto A^. This construction was introduced in [22]. The interest 
of this construction is that it provides a Riemannian metric on the image manifold A^ 
which incorporates the group actions of G. The evolution equations were derived and 
studied in ^24j in the case liu^n^p) = WuW^ + ||z^||^, for a given norm, || . ||g, on q and a 
prescribed Riemannian structure on the manifold A^. In |17| the metamorphosis approach 
to image matching was formulated in terms of Hamilton's principle in the Euler-Poincare 
variational framework |T2] to derive the evolution equations. This formulation provides 
an interesting contrast between the variational formulations of optimal control problems 
and evolutionary equations. 




Li9,9,V,v) = L{gg \gv,9V) =■ Ku,n,u) = ||(n,z/)||^ 
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Semidirect product Lie groups Assume that iV is a Lie group and that for all 
g E G, the action of on is a group homomorphism. That is, for all n,h G N, 
g{nh) = {gn){gn). For example, N can be a vector space and the action of G can be 
hnear. Consider the semidirect product G(§)N with {g,n){g,n) = {gg, {gn)n) and build on 
G(§)N a right-invariant metric constrained by its value || . ||(idG,id]v) at the identity. Then 
optimizing the geodesic energy in G(§)N between (idc^o) and {gi,ni) with fixed uq and 
rii and free gi yields a particular case of metamorphosis. 

The Euler-Poincare formulation of metamorphosis on Diff(s)G was also considered in 
[T7] . The Euler-Poincare equations for a Lagrangian /(m, z/) in which the variable n is 
absent are found for G(S)N = Diff(s)jF to produce the CH2 and MCH2 systems for the 
corresponding norms, 

lcH2{u, v) = \\u\\%i + ||z/||^2 and Imch2{u, v) = \\u\\%i + . 

These are precisely the norms discussed in the present paper. Thus, the theory of meta- 
morphosis in imaging science summons the singular solutions for MCH2, when the smooth- 
ing kernels Ki and K2 are chosen to be Helmholtz inversions. Future developments will 
determine the utility of these singular solutions in imaging science. 

Actually, the problem of matching shapes in imaging science is an optimal control 
problem, rather than the initial value problem (IVP) discussed here. However, as we have 
seen, the IVP for MCH2 is illuminating and interesting in its own regard. Moreover, 
considerations of the IVP help in the interpretation of the solutions of the optimal control 
problem for image matching as a flow of information parcels (landmarks) from one image 
to another. The particle-hke solution behavior for the IVP discussed here is expected to 
hold for any symmetric, confined, translation-invariant choices of the smoothing kernels 
Ki and K2. 
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